      subroutine rgefa(a,lda,n,ipvt,info)
      integer lda,n,ipvt(n),info
      real a(lda,n)
c
c     dgefa factors a real matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       real(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c
c     internal variables
c
      real t
      integer idamax,j,k,kp1,l,nm1
c
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0e0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0e0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0e0) info = n
      return
      end
      subroutine  dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real da,dx(n)
      integer i,incx,m,mp1,n,nincx
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      integer function idamax(n,dx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      real dx(n),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if( n .lt. 1 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      dmax = abs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(abs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = abs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 dmax = abs(dx(1))
      do 30 i = 2,n
         if(abs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = abs(dx(i))
   30 continue
      return
      end
c$$$      subroutine daxpy(n,da,dx,incx,dy,incy)
c$$$c
c$$$c     constant times a vector plus a vector.
c$$$c     uses unrolled loops for increments equal to one.
c$$$c     jack dongarra, linpack, 3/11/78.
c$$$c
c$$$      real dx(n),dy(n),da
c$$$      integer i,incx,incy,ix,iy,m,mp1,n
c$$$c
c$$$      if(n.le.0)return
c$$$      if (da .eq. 0.0e0) return
c$$$      if(incx.eq.1.and.incy.eq.1)go to 20
c$$$c
c$$$c        code for unequal increments or equal increments
c$$$c          not equal to 1
c$$$c
c$$$      ix = 1
c$$$      iy = 1
c$$$      if(incx.lt.0)ix = (-n+1)*incx + 1
c$$$      if(incy.lt.0)iy = (-n+1)*incy + 1
c$$$      do 10 i = 1,n
c$$$        dy(iy) = dy(iy) + da*dx(ix)
c$$$        ix = ix + incx
c$$$        iy = iy + incy
c$$$   10 continue
c$$$      return
c$$$c
c$$$c        code for both increments equal to 1
c$$$c
c$$$c
c$$$c        clean-up loop
c$$$c
c$$$   20 m = mod(n,4)
c$$$      if( m .eq. 0 ) go to 40
c$$$      do 30 i = 1,m
c$$$        dy(i) = dy(i) + da*dx(i)
c$$$   30 continue
c$$$      if( n .lt. 4 ) return
c$$$   40 mp1 = m + 1
c$$$      do 50 i = mp1,n,4
c$$$        dy(i) = dy(i) + da*dx(i)
c$$$        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
c$$$        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
c$$$        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
c$$$   50 continue
c$$$      return
c$$$      end
